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SUMMARY 

The  computer  program  ZMODE  computes  the  eigenfunctions  and  dispersion 
relations  for  internal  wave  oscillations  in  a density-stratified  ocea., 
thermocline.  The  computation  is  rapid  and  reliable,  owing  to  a 
procedure  which  obtains  very  accurate  trial  eigenvalues  by  exploiting 
certain  analytic  properties  of  the  differential  eigenfunction  equation. 

This  document  describes  in  detail  the  analytic  procedures  and  program 
structure,  and  provides  complete  operating  instructions. 
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Because  the  derivatives  are  obtained  analytically,  via  properties  of 
equation  (1),  rather  than  numerically,  the  two  tables  can  be  used  to 
construct  highly  accurate  interpolated  dispersion  data.  For  example, 
relative  interpolation  accuracies  better  than  10"6  have  been  obtained 
along  dispersion  curves  represented  by  thirty  points. 

A particularly  efficient  and  reliable  procedure  is  used  for  the 
eigenvalue  search,  in  which  the  derivatives  dyadic  are  used  to  extrapolate 
successive  trial  values  along  each  eigenlocus.  Converged  eigenvalues  are 
consequently  obtained  in  just  over  two  iterations  on  the  average,  and 
the  computation  time  per  eigenvalue  is  (CDC  7600) 

( eqn  (1)  integration  steps/100)  x 0.5  x 10~2  sec. 

ZMODE  consists  of  two  modules:  TCLINE,  which  converts  input 

temperature-versus-depth  data  to  the  stability  profile  N2(z),  and  ZMODES, 
which  performs  the  eigenfunction  calculations.  The  dummy  module  UTIL  is 
included  for  the  convenience  of  the  user,  who  can  replace  it  with  sub- 
routines performing  various  calculations  on  the  quantities  provided. 

Each  module  provides  its  own  input  and  output  functions,  and  operates 
essentially  as  an  independent  program.  Input,  output  and  operating  parameter 
specifications  are  communicated  directly  to  the  module  named  on  the  control 
card.  While  this  level  of  organization  may  seem  a little  elaborate,  it 
has  been  imposed  with  the  idea  that  ZMODE  will  ultimately  be  the  nucleus 
of  much  larger  internal  wave  codes.  It  therefore  seemed  desirable  to 
establish  a modular  organization  and  uniform  control  procedure  so  that  new 
program  modules  can  be  added  and  used  as  simply  as  possible.  For  this 
reason  also  the  various  kinds  of  numerical  and  graphical  output  are  optional, 
and  must  be  specifically  requested  by  appropriate  control  cards. 


-3- 


2.  Sample  Calculation 

The  types  of  numerical  and  graphical  output  available  from  ZMODE 
are  illustrated  in  the  pages  following  by  reproductions  from  a test 
calculation.  The  shallow-water  geometry  and  thermal  scale  of  this 
particular  case  are  characteristic  of  the  NUC  to  er  site  near 
San  Diego,  while  the  particular  double  thermocline  profile  used 
is  a fictitious  one  designed  to  exhioit  intermode  resonances,  as 
will  be  seen. 

Figure  1 shows  the  numerical  output  of  the  interpolated  tempera- 
ture and  stability  profiles,  and  Figures  2 and  3 show  the  graphical 
counterparts. 

The  dispersion  calculations  for  this  case  encompassed  modes 
1-5,  each  on  a wavenumber  range  of  0 to  280  cycles/km  in  41  increments. 
Execution  time  on  the  CDC  7600  was  1.25  seconcs  for  the  dispersion 
calculations  themselves  and  1.7  seconds  overall.  The  detailed 
dispersion  printout  includes  frequency,  phase  speed  and  group 
speed  for  each  mode,  as  shown  in  Figure  4. 

Graphical  display  of  the  dispersion  curves  is  available  as  in 
Figure  5.  Note  the  "Eckert  Resonance"  between  modes  2 and  3 at 
145  cycles/km,  where  nearly  independent  oscillations  on  the  two 
stable  layers  become  strongly  coupled  because  of  ar  accidental 
coincidence  in  frequency.  The  relative  frequency  difference  at 
closest  approach  is  0.013,  which  is  not  resolvable  on  the  printer 
plot. 

Printer  plots  of  the  eigenfunction  profiles,  as  shown  in 
Figures  6-8,  are  useful  in  interpreting  thermocline  behavior. 

The  six  curves  shown  for  each  mode  correspond  to  six  equally-spaced 
wavenumbers  from  zero  on  the  left  to  maximum  on  the  right. 


Note  how  at  high  wavenumber  the  oscillations  are  confined  principally 
to  one  of  the  two  stable  layers.  Note  also  the  sudden  exchange 
occurring  in  principal  layers  between  modes  2 and  3 at  the  resonance. 
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***  THEHHOCLlNt  PROFILE  *** 


DEPTH,  M 

TtMP,  C 

NCZ),  CY/HR 

I 

-0,00 

20,00 

0.00 

-.40 

20,01 

0,00 

-,80 

20,01 

3,80 

• 1,20 

19,99 

6,99 

-1,60 

19,96 

8,71 

• 2,00 

19,91 

9,60 

• 2,«0 

19,06 

9,91 

-2,80 

19,01 

10,15 

• 3,20 

19, 75 

11,37 

• 3,60 

19,67 

13,34 

-a, oo 

19,56 

15,59 

• 4,40 

19,41 

17.24 

• 4,80 

19,24 

18,32 

• 5,20 

19.05 

19,28 

• 5,60 

18,82 

21,90 

• 6,00 

18,50 

25,81 

• 6,40 

18,07 

27,58 

• 6,80 

17,65 

25,08 

• 7,20 

17,51 

22.07 

• 7,60 

17,03 

20,23 

• 8,00 

16,79 

18,48 

• 8,40 

16,60 

16,23 

• 8,80 

16,95 

13,23 

• 9,20 

16,56 

9,74 

• 9 , <>0 

16,52 

6,32 

•10,00 

16,30 

3,08 

•10,40 

16,30 

1,72 

•10,80 

16,29 

3,65 

•11,20 

16,28 

6,14 

•11,60 

16,24 

8,71 

•12,00 

16,17 

1 1,28 

•12,40 

16,07 

13,59 

■12,80 

15,92 

15,20 

•13,20 

15,75 

16,31 

• 13  ,t>0 

15,55 

16,90 

•14,00 

15,34 

16,76 

•14,40 

15,15 

15,93 

•14,80 

14,97 

14,53 

•15,20 

14,03 

12,99 

•15,60 

14,71 

11.27 

•16,00 

14,63 

9,28 

•16,40 

14.58 

7,51 

•16,80 

14,54 

6,46 

•17,20 

14,51 

6,46 

•17,60 

14,47 

7,17 

•18,00 

14,43 

7,68 

•18,40 

14,38 

7,99 

•10,80 

14,33 

8,13 
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Figure  1 

Numerical  Thermocline  Profile 


Figure  2.  Graphical  Temperature  Profile 
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**«  MODE  1 DISPERSION  RELATIONS  *** 


K,  CY/KH 

h,  1/SEC 

CY/HR 

C,  M/SEC 

CG,  6/SEC 

0,00000 

0,00000 

0,0 JOOO 

,16912 

,16912 

7,00000 

.00721 

4,12816 

.16382 

i 15375 

14,00000 

.01325 

7.594  19 

,15066 

,12014 

21.00000 

.U  1 780 

10.19563 

, 1 3487 

,08753 

28,00000 

.02109 

12,08173 

, 1 1906 

,06368 

35,00000 

.02352 

13.47359 

,10  f 3 

,04792 

42,00000 

.02538 

14,54258 

,09610 

,03763 

49,00000 

,02687 

15,39780 

.08729 

,03068 

56,00000 

.02811 

16,10555 

,07989 

,02575 

63,00000 

,02916 

16.70603 

.07366 

,02208 

70,00000 

,03006 

17,22531 

,06035 

,01925 

77.0COOO 

.03086 

17,60097 

,06378 

,01700 

84,00000 

,03157 

10,08556 

,05901 

,01517 

91,00000 

,03220 

10,44831 

.05631 

,01366 

98,00000 

.03277 

18,77626 

,05322 

,01240 

105,00000 

,03329 

19,07484 

,05046 

, C 1133 

112,00000 

.03377 

19.34837 

,04799 

,01040 

119,00000 

.03421 

19,60027 

,04575 

,00961 

126,00000 

.03462 

19,83335 

,04372 

,00091 

133,00000 

,03494 

20,0499  1 

,04188 

,00029 

140,00000 

,03535 

20.25187 

,04018 

,00775 

147,00000 

,03568 

20,44082 

,03863 

,00726 

154,00000 

.03599 

20,61815 

,03719 

,00602 

161,00000 

, 03628 

20.78502 

, 0 350C 

,00643 

168,00000 

,03655 

20,94245 

.03463 

,00607 

175.00000 

,03681 

21,09130 

,03340 

,00575 

182,00000 

,03706 

21,23236 

,03241 

,00545 

189,00000 

,03729 

2 1 ,36629 

,03140 

,00518 

19b, 00000 

.03751 

21 ,49369 

, 03046 

,00493 

203,00000 

, 03773 

21 ,61500 

,02958 

,00470 

210,00000 

,03793 

21,73094 

,02874 

,00449 

217,00000 

,03012 

21,84168 

,02796 

,00430 

224,00000 

,03631 

21,94769 

, 02722 

,00412 

231,00000 

,03048 

22,04929 

,02651 

,00395 

238,00000 

,03865 

22,14679 

,02585 

,00379 

245,00000 

,03602 

22,24047 

,02522 

,00365 

252,00000 

,03897 

22,33058 

,02461 

,00351 

259,00000 

,0391  S 

22,41735 

,02404 

,00338 

266,00000 

,03927 

2?,50097 

,02350 

,00326 

273,00000 

,03991 

22,58105 

,02298 

,00315 

280,00000 

,03955 

22,65955 

,02248 

,00304 

Figure  4.  Numerical  Dispersion  Relations 
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Figure  8 

Mode  3 Eigenfunction  Graphical  Profile 
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3-  Operating  Instructions 


All  program  functions  are  governed  by  control  cards  in  a standard  format 
of  three  ten-character  hollerith  instruction  words  (left-justified)  followed 
by  up  to  five  ten-character  numerical  parameters: 


INPUT 


TCLINE 
| ZMODES I 


PRINT 


(Option 

descriptor) 


I < 


<_  5 parameters 


GRAPH 


The  first  word  describes  the  function  to  be  performed  and  the  second  describes 
the  module  to  which  the  instruction  is  addressed.  The  instructions  are  self- 
explanatory;  the  distinction  between  INPUT  and  SET  is  that  INPUT  refers  to 
tabular  data  and  SET  refers  to  adjustable  program  parameters.  Cards 
refering  to  different  modules  can  be  intermixed,  the  only  constraints  being 
the  obvious  ones  that  INPUT  and  SET  must  precede  the  first  EXEC  for  a 
given  module,  and  PRINT  and  GRAPH  must  follow.  EXEC  TCLINE  must  precede 
EXEC  ZMODES.  Any  number  of  control  cards  can  be  used.  For  example, 
repeated  calculations  with  varying  thermocline  data,  wavenumber  increments, 
or  wavenumber  range  can  be  requested,  and  duplicate  output  for  a given 
calculation  can  be  obtained  with  multiple  output  cards. 


Cards  with  arbitrary  alphanumeric  content  placed  before  the  first 
control  card  will  be  duplicated  on  the  output  title  page  and  otherwise 
ignored;  these  can  be  used  as  desired  for  descriptive  titles.  Execution 
is  terminated  by  a card  with  a left-justified  STOP.  The  various  control 
parameters  are  described  below. 


I. 


I 


I 
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INPUT  TCLINE:  (NDATA) 

Format  ( 110  ) 

This  instruction  is  followed  by  cards  containing  the  temperature- 
versus-depth  data  in  order  of  increasing  depth, 

Zi > i = 1,  2,  ....  NDATA  units:  Meters , Cel  ci  us 

in  the  format  (2F10.J.  Z]  must  be  zero.  Z (NDATA)  must  be  greater  than 
or  equal  to  the  parameter  ZA,  described  in  the  following.  NDATA  must  be 
in  the  range  (3,100). 

SET  TCLINE:  (ZA,  ZB,  N ) 

Format  (2F10._,  110) 

ZA  is  the  depth  (M)  of  the  themocline  bottom,  below  which  it  is  assumed 
that  density  stratification  vanishes;  ZA  must  be  less  than  or  equal  to 
the  depth  of  the  last  input  data  point.  ZB  (M)  is  the  depth  of  the  ocean 
bottom,  greater  than  or  equal  to  ZA.  The  parameter  N defines  the  number 
of  integration  steps  between  zero  and  ZA  used  in  t^e  mode  calculations 
and  is  specified  here  so  that  the  interpolated  stability-frequency  table 
created  by  TCLINE  will  have  the  right-sized  increments.  A useful  rule 
of  thumb  is  that  N should  exceed  ten  times  the  highest  mode  order  to  be 
calculated;  maximum  permissible  value  is  200. 

EXEC  TCLINE:  (No  parameters) 

PRINT  TCLINE:  (ISKIP) 

Format  ( 110  ) 

Up  to  399  values  (i.e.,  2N-1)  of  interpolated  temperature  and  stability 
frequency  can  be  printed.  For  shortened  printout,  only  every  ISKIPth_ 
value  is  printed.  When  ISKIP  is  blank  the  default  value  4 is  used.  See 
Figure  1 , p. 5. 
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GRAPH  TCLINE:  (T1  , 12 , WSCALE) 

Format  ( 3 FI 0 . ) 

This  produces  a printer  plot  of  T in  the  interval  (T2,  T1 ) versus  Z in  the 
interval  (0,ZA),  and  a plot  of  N(Z)  in  the  interval  (0,  WSCALE)  versus  Z. 

Units  are  meters,  C,  and  cy/hr.  See  the  example  in  Figures  2 and  3,  pp.  6&7. 
If  the  proper  value  of  WSCALE  is  not  known,  a rough  value  can  be  estimated 
and  several  values  tried  at  once  with  repeated  GRAPH  TCLINE  cards. 

INPUT  ZMODES : {Not  Used) 

SET  ZMODES:  ( KMAX , NK,  ERR) 

Format  { FI 0 110,  E10.0) 

This  prepares  ZMODES  to  find  eigenfunctions  and  eigenvalues  at  NK  ( <51 ) 

equally-spaced  values  of  wavenumber  k in  the  interval  0 <_  k £ KMAX. 

The  units  of  KMAX  are  cycles/km,  converted  internally  to  radians/m. 

ERR  specifies  the  minimum  relative  precision  to  which  the  iterated 

eigenvalue  search  must  converge.  The  default  value  is  10"8,  but  up  to 
*“13 

10  is  usable  on  a 60  bit  machine. 

EXEC  ZMODES:  {Ml , M2,  Z1 , Z2) 

Format  (2110,  2F10.J 

The  above  instructs  ZMODES  to  find  eigenfunctions  and  eigenva^es  for 
modes  Ml  through  M2  <_  20.  As  calculations  are  completed  for  each  mode 
a short  line  of  diagnostic  output  is  produced  listing  the  required  computation 
time  and  average  number  of  search  iterations  per  eigenvalue.  See  Figure  9. 

Two  depths,  Z1  and  Z2,  can  be  specified  in  the  interval  (0,ZA),  at  which 
eigenfunction  values  and  Z-deri vatives  will  be  tabulated  along  with  dis- 
persion data  during  execution  of  ZMODES. 


***  tIGENVL  COMPUTAT ION  TIMES  *** 


* RESOLUTION  41  VALUES/MODE  * 

i 


MOCE 

AVG  ITERATIONS 

TIME 

i 

i 

2,71 

,2350 

2 

2,83 

,2480 

3 

3.20 

,2760 

41 

2,66 

,2340 

5 

3,12 

,2690 
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PRINT  ZMOOES: 
Format 


(Ml  , M2) 
( 2110  ) 


Provides  printed  tabular  output,  versus  wavenumber,  of  frequency,  phase 
speed,  and  group  speed,  for  modes  Ml  through  M2.  See  Figure  4,  p.8. 


GRAPH  ZMODES 


Format 


j Ml,  M2,  WSCALE,  KSrALE  I 

FUNCTIONS I )M1 , M2,  MULT,  FSCALE,  ZSCALE  ) 


( I10’  2F10.J 


! 


With  the  blank  option  this  instruction  plots  the  frequency  versus  wave- 
number  dispersion  curves  for  modes  Ml  through  M2  as  shown  in  Figure  5, 

P.  9.  Plot  ranges  are  (0,  KSCALE)  cy/km  in  wavenumber  and  (0,  WSCALE) 
cy/hr  in  frequency. 

The  FUNCTIONS  option  plots  the  eigenfunctions  themselves  side  by  side 
for  MULT  (>_  2)  equally  spaced  wavenumbers  between  0 and  KMAX , one  mode 
per  figure  for  modes  Ml  through  M2.  Figures  6-8,  PP. 10-12,  show  examples 
for  modes  1-3  with  MULT=6. 

The  functions  a. e automatically  normalized  for  plotting  and  the  dimension- 
less number  FSCALE  is  available  for  further  scaling.  A fair  rule  of 
thumb  for  FSCALE  is 

FSCALE  MULT  X M2 


but  this  can  vary  with  thermocline  character. 
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Stable  eigenvalue  searches  along  the  eigemoci  require  that  the 
extrapolated  eigenvalue  predictions  be  accurate  relative  to  the  inter- 
mode eigenvalue  spacing.  At  large  wavenumbers  this  stability  can  be 
jeopardized  by  twc  effects:  long  exponential  scales  in  the  eigenfunction 

integration  introduce  numerical  noise  into  the  calculation  of  eigenlocus 
slope  and  degrade  the  extrapolation  accuracy;  so-called  "Eckart  resonances" 
between  isolated  stable  layers  become  severe  and  intermode  eigenvalue 
spacing  can  become  exponentially  small.  In  this  context  wavenumbers  are 
"large"  or  "small"  in  the  comparison 

KMAX  * 1300  . 

fi  ft 

Convergence  precision  (ERR)  of  10  -10  will  probably  be  adequate  for 

small  wavenumbers,  and  10“^  for  the  larger  wavenumbers  of  practical 
interest.  In  limited  tests  and  ERR  setting  of  10“13  removed  all  evidence 
of  numerical  noise  out  to  KMAX  • ZA  = 10,000. 

Since  the  error  in  the  quadratic  eigenvalue  extrapolation  scale:  as 
_ 3 

NK  , the  ability  c*  the  extrapolation  to  get  through  an  Eckart  resonance 
can  be  improved  sharply  by  an  increase  in  NK. 


An  accidental  convergence  to  the  wrong  mode  will  result  in  program 
termination  with  the  error  message  "EIGENVL  failed  to  converge  at 
mode , wavenumber . " 


r ; , 

La 

#v  .< 


‘1 


The  average  number  of  iterations  per  eigenvalue  necessary  to  converge 
to  ERR,  printed  for  each  mode  during  ZM0DES  execution,  is  an  indication 
of  the  quality  of  :trapolations  being  obtained;  less  accurate  trial 
eigenvalues  require  more  iterations.  The  number  should  be  just  above 
two  for  small  KMAX,  large  NK,  and  ERR  ^ 10  Values  around  three  are 
to  be  expected  for  ERR  < 10  Since  total  execution  time  is  proportional 


J 


I 
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to  NK  times  the  average  number  of  iterations,  the  smallest  NK  and  largest 
ERR  consistent  with  reliable  computation  should  be  aimed  for  when 
computing  expense  is  a factor. 
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5.  Utility  Subroutine 


Various  quantities  calculated  by  ZMODE  are  available  to  the  user  in 
the  dummy  subroutine  UTIL: 


r 


I 

i 


I 

! 

f 


f. 


Quantity 

Symbol 

FORTRAN 

Wavenumber 

k 

TK  (IK) 

Inverse  phase  speed 

Y (k) 
Imv 

TGAMMA  (IK,  MODE) 

Associated  derivative 

Vdk 

TDGAMMA  (IK,  MODE) 

Mode  amplitude 

»m<zj  > k> 

FI  (IK,  MODE,  JZ) 

Associated  derivative 

3<j)  / 3Z 

FIZ  (IK,  MODE,  JZ) 

Mode  sample  depths 

V1*2 

Zl,  12 

Table  lengths 

(IK  ■ 1,  NK) 

NK 

From  the  first  three  tables 

one  can  obtain  a cubic 

interpolation  of 

7m  for  any  k.  Related  dispersion  quantities  can  then  be  computed  via 


c = y"1  (phase  speed) 

m 'm 

w = ky'1  (frequency) 


= daim/dk 


= c (1  - ui  dy  /dk ) (group  speed); 
m mm 

the  units  are  (meters,  seconds,  radians). 


I 
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The  mode  amplitude  and  derivative  data  are  retained  at  only  two  depths, 
Z1  and  Z2,  as  specified  cn  the  EXEC  ZMODES  card.  These  are  intended  fov 
use  in  constructing  two-point  Green's  functions  needed  in  studies  of 
internal  wave  generation.  Another  possible  application  is  the  pre- 
diction of  surface  current  components  u^  (0,  k)  associated  with  measured 
isotherm  displacements  (z-j)  at  a given  depth  z^: 

f -i  . d*m(0)l 


Ym  ^zl^ 


(5-1) 


* 


6.  Methods  of  Calculation 


The  thermocline  is  numerically  represented  by  a table  of  stability 

frequency  N (z)  for  2n-l  equally-spaced  depths  spanning  the  interval 

(-za,  0);  outside  this  interval,  that  is,  in  the  remaining  interval 
a 2 

(-zb,  -za)  to  the  ocean  bottom,  N (z)  is  assumed  to  vanish.  The  numbers 
za,  zb  and  n are  input  parameters.  The  table  is  obtained  from  an  input 
table  of  temperature  data,  T(z.j),  by  a third-order  spline  interpolation 
and  differentiation,  via  the  formula 

N2(z)  = - ct(T)  g |jl  (6-1) 

in  which  a(T)  is  a linear  representation  of  the  coefficient  of  thermal 
expansion  of  seawater  at  35%  salinity.  The  use  of  the  third-order  spline 
yields  a smooth  stability  profile  with  a continuous  derivative. 

The  eigenvalue  equation 

+ (Y2  N2  - k2)  f =*  0 , (6-2) 

dz^ 

here  written  in  terms  of  wavenumber  k and  inverse  phase  speed  y as 
parameters,  has  solutions  that  satisfy  the  boundary  conditions 

f(-zb)  = f(0)  = 0 (6-3) 

on  an  infinite  number  of  real  loci  Ym(k)-  At  each  k,  the  eigensolutions 
obey  the  orthogonality  rule 
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/ f N2  f dz  = u <5  . (6-4) 

J7  m n mn  1 ' 

"Id 

The  numerical  search  for  the  eigenfunctions  f , eigenvalues  y , and 
normalizing  factors  um  proceeds  as  follows.  For  a given  value  of  k 
and  a trial  value  of  y assumed  close  to  the  rr>^-  eigenlocus,  the  trial 
function 


f = f(z,  y , k),  f(-2J  = 0 

a 

is  numerically  integrated  in  n steps  from  -z  to  0,  the  appropriate 
starting  values  of  fm(-za)  and  f^(-za)  having  been  found  from  the 
known  analytic  solution 


f = const,  x sinh  k (z  + zb)  (6-5) 

in  the  region  (-^,  -z^)  where  N = 0.  The  value  of  f at  z = 0, 

w(y,  k)  = f(0,  y,  k)  (6-6) 

is,  in  general,  not  zero  and  vanishes  only  wheny=  y (k).  The  quantity  w 
is  a continuous,  differentiable  function  of  y and  k,  and  its  properties 
are  key  to  the  eigenvalue  search  procedure.  If  we  define 

g(z,  Y , k)  = (z,  y , k)  (6-7) 

3y  ^ 

and  formally  differentiate  the  eigenvalue  equation  with  respect  to  y , 
we  see  that  g satisfies  the  inhomogeneous  equation 


|W| 


m 


. CSr. 
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+ (y2  N2  - k2)g  = - N‘f,  (6-8) 

dzd 


g(-2a)  = g'C-Zg)  = 0. 

This  companion  equation  is  simultaneously  integrated  with  (6-2),  yielding 


1 3w 
2Y  3Y 


= g(o , y » k)  , 


(6-9) 


which  can  be  used  immediately  to  generate  an  improved  estimate  of  y, 


y ■+  Y 


(f)'1  - • 


(6-10) 


A few  such  iterations  drive  w quickly  to  zero,  within  the  limit 
of  computational  precision.  The  numerical  integration  of  the  similar 
equations  (6-2,  6-8)  uses  constant  step  size  and  identical  algorithms. 
Since  the  resulting  finite  difference  equation  for  g is,  in  fact,  the 
formal  derivative  of  the  finite  difference  equation  for  f,  the  con- 
vergence is  very  smooth,  and  the  correction  in  (6-10)  will  approach 
within  a few  bits  of  zero. 


The  value  of  3w/3y  has,  upon  convergence  to  the  eigenvalue,  the 
following  significance.  Multiplying  (6-2)  by  g , (6-8)  by  fm,  and 
subtracting  , we  have 


d 

37 


dg\ 

rm  37~y 


■ - N2  fjji 

m 


(6-11) 


l 
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i n tegra ting  this  expression  and  remembering  that  fm  vanishes  at  -z 


and  0,  while  g vanishes  at  -z^,  we  find 


— ) f ' (0)  = -f°  N2  f2  dz  = . 

3Y / m v J _ m m 


(6-11  ) 


Thus,  the  calculated  quantity  3w/3y  furnishes  the  normalization  constant 


y . An  analogous  argument  shows  that 


J_  (M.)  fy(o)  = v 
2k  V ak  / V ; r 


(6-12) 


where  v is  the  simple  norm 
m r 


vm  = / r dz 

m */  7 m 


(6-13) 


Upon  completion  of  each  iteration,  this  last  quantity  is  calculated  by 
direct  integration  of  the  squared  eigenfunction.  The  slope  of  the  eigen- 
locus,  dYm/dk,  defined  by 


K I 


is  then  obtained  as 


W dYm  + dk  = °’ 


'm  _ 3w/3k 


3w/,lrm  Vm 


(6-14) 


Formally,  the  identity  above  is  a differential  equation  for  the  entire  m- 
eigerlocus.  Numerically,  it  is  used  to  extrapolate  accurate  trial  eigen- 
values at  successive  points  along  the  locus, 


Ym  (k+Ak)  = Ym(k)  +[?  IT  (k)  " 2 alf  (k"Ak)J  Ak 

+ O(Ak)3  . ( 

Starting  trial  eigenvalues  for  all  the  loci  at  k=0  are  provided  by  the 
WKB  approximation  to  (6-2). 

Ym  (0)  J°  N (z)  dz  ^ (m  - n ; ( 

"za 

the  first  extrapolation  from  k=0  is  modified  from  (6-14)  to  use  the 
second  derivative 

dk^  Ymym 


obtained  from  (6-13)  by  L'Hopital's  rule. 


i 


7.  Program  Structure  and  Nomenclature 

The  subroutine  calling  sequence  and  data  flow  are  depicted  in 
Figure  10.  The  output  data  tables  comprising  COMMON/FLD/  have  already 
been  described  in  section  5.  Other  variables  and  working  files  are: 


Description 

| I 

Symbol 

FORTRAN 

Integration  step  size 

dz 

DZ 

Number  of  steps 

n 

N 

Stability  profile 

N2(z) 

QN(X), 

1=1,  N2 

2n-l 

N2 

Trial  eigenfunction 

f(z, Y , k) 

F ( I ) 

1=1,  N 

I I 

( 3f/ 9z)  • dz 

FZ(I) 

Associated  function 

g(z,  y , k)/dz2 

6(1) 

1=1,  N 

(3g/3Z)/dZ 

GZ(I) 

Wavenumber 

k 

K (real) 

Increment 

Ak 

DK 

Trial  phase  slowness 

Y 

GAMMA 

Temperature  profile 

T(z) 

TDATA (1)1 

1=1 

ZDATA(I)  j 

Mode  number 

m 

MODE 

i 


8-  Description  of  Routines 


a.  ZMODE 


ZMODE  is  the  driver,  whose  only  function  is  to  pass  control -card  images 
to  the  module  designated:  TCLINE,  ZMODES,  or  UTIL. 


b.  TCLINE 


TCLINE  reads  the  temperature-versus-depth  data,  calls  SPLINES  to  compute 
the  cubic  spline  coefficents,  then  uses  SPLINE  (ENTRY  SPLINE2)  to  obtain 
interpolated  values  of  dT/dz  at  2n-l  equally-spaced  points  between  -z 

d 

and  0.  The  stability  profile  is  computed  and  tabulated  as 


N2(z)  = -g(Cl+c2T)  £ 


(8-1) 


in  which  c-|+c2T  is  a linear  fit  to  the  coefficent  of  thermal  expansion  of 
water  at  35"  salinity.  TCLINE  also  computes  the  quantity. 


N(z}dz 


(8-2) 


for  later  use  by  ZMODES  in  obtaining  starting  eigenvalues  via  the  WKB 
approximation  (see  egn.  8-16). 


c.  ZMODES 


ZMODES  computes  the  dispersion  tables  y (k)  and  d ym/d k for  all  modes  in 


in  the  interval  (Ml,  M2)  designated  on  the  execution  card. 


For  each  mode,  ZMODES  computes  a starting  trial  eigenvalue  for  k = 0 via 


the  WKB  approximation,  calls  EIGENVL  to  compute  ym  and  dym/d k to  the 


specified  precision,  and  thereafter  repeatedly  increments  k,  calls  EIGENVL, 
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and  taoulates  y , dy  /dk  up  to  the  specified  wavenumber  limit.  Successive 
m m r 

trial  eigenvalues  are  computed  according  to  the  scheme  described  in  section 
6 and  passed  to  EIGENVL  with  each  call.  ZMODE  will  terminate  execution 
with  an  error  message  if  the  diagnostic  integers  IERR  and  MERR  returned  by 
EIGENVL  are  not  both  zero  (sto  below). 


For  the  printer  plot  of  the  dispersion  curves  “m(k)  an  auxiliary 
200-point  table  is  computed  in  convenient  cy/hr  units  from  cubic  inter- 
polations of  y (k)  provided  by  FINT.  For  the  printer  plots  of  the  un- 
normalized eigenfunctions  a scale  factor  derived  from 


v 

m 


/<> 


is  passed  to  the  plotting  subroutine  GRAPH. 


d.  EIGENVL 


This  subroutine  provides  the  iteration  logic  for  the  convergent  eigenvalue 
search  at  a given  mode  and  wavenumber,  and  computes  the  eigenfunction  norm 
parameters  used  by  ZMODES  for  trial  eigenvalue  extrapolation.  The  argument 
list  is 


GAMMA  Inverse  phase  speed  y; 

trial  value  on  input, 
converged  value  on  return 

K Wavenumber  k 

NORM  Norm  v 

ONORM  Weighted  norm  u 


MODE 


Desired  mode  number 


MERR 


Mode  error  flag 


IERR  Iteration  failure  flag 

IT  Number  of  iterations  required 

EIGENVL  repeatedly  calls  EIGENFN  to  perform  integrations  of  the  eigenvalue 
equation  with  the  current  value  of  y,  each  time  improving  y by  Newton's 
method 


until  w is  sufficiently  close  to  zero  that 

I Ay/y|  < ERR. 

Accidental  convergence  to  an  adjacent  mode  is  prevented  by  a check  of  the 
sign  of  dw/dy,  which  alternates  as  a function  of  mode  number.  If  the  sign 
is  wrong,  or  if  for  some  reason  convergence  has  not  been  obtained  after 
ten  iterations,  IERR  is  changed  from  0 to  1 . 

The  number  of  zero  crossings,  M,  occuring  in  f(z)  (excepting  z-o)  is 
counted  and  the  number 

MERR  = M - MODE  + 1 

formed;  if  either  IERR^O  or  MERR^O  further  computation  is  bypassed. 

After  normal  convergence  EIGENVL  makes  a linear  correction  to  the  eigenfunction 
using  the  last  computed  value  of  Ay, 


f (z)  ► f(z)  - 2yAyg(z). 


(3-3) 
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The  reason  for  this  correction  is  that  a small  relative  error  in 

eigenvalue,  Ay/y,  can  be  exponentially  amplified  by  the  eigenvalue 

2 2 

equation  in  regions  where  > N (z),  leading  to  errors  in  f of 
order 


£f  ,,  iL  ekza 

f y ’ 


so  that,  for  example,  at  kz  ^ 30  the  error  can  grow  to  order  unity 

-13  d 

even  for  Ay/y  = 10  . Fortunately,  the  error  is  predictable  and 

removable  by  (8.3),  to  order 


M2  Mi 


(8-4) 


The  norm  quantity  ym  is  then  computed  via 


- - X rm  to). 

2 'I'm  s'%  m 


and  the  quantity  by  integration, 


fm2  dz' 


analytically  in  the  region  (-zb>  -za),  and  numerically  according  to  the 
end-corrected  trapezoidal  rule  in  the  region  (-z  ,o),  with  a small  predicted 
error  of  order 

2 (4)  4 

za(n  (dz) 


e.  EIGENFN 


EIGENFN  integrates  f and  the  associated  function  g on  the  interval 


L 


(-za,o),  using  a fourth  - order  algorithm  furnished  by  G.  Peebles  of  RDA: 
yi+l!  - y,  + 7 Azy!  ♦ 1 AZ2yV 

* y"  (y,-+s) 


^i+>!  + 24AZ  (y'i+'i 


(M 


y,  + Azyl  + 1 az2  (yV  + 2yf”+lj) 


1+1  1 y 


yi+i  =yi  +r42  (y'i + 4yiH  + yi'+i) 


A constant  step  size  az  is  used  to  allow  convenient  non-dimensional ization, 
with  the  consequent  savings  in  arithmetic. 

The  arguments 

W = -f(0) 


m 


WG  - -g(0) 


are  set  prior  to  return. 


r 

J- 


i 

t 


,c 

i . 


f. 
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F I NT 


The  function  subroutine  PINT  performs  a cubic  interpolation  of  a tabulated 
function  with  the  aid  of  corresponding  tabulated  derivative.  The  arguments 
are 

Z Interpolation  point 


i 


r 


ZT  Array  containing  values  of  independent  variable  Z 

FT  Array  containing  values  of  tabulated  function 

FTZ  Array  containing  values  of  tabulated  derivative 

N Dimension  of  arrays 

I Pointer  index 

The  subroutine  logic  finds  I such  that  ZT ( I ) <_  Z < ZT ( I +1 ) prior  to 
interpolation,  and  the  value  of  I can  be  held  for  subsequent  calls  to 
minimize  table  searches. 


g.  SPLINES 

SPLINES  computes  third-order  spline  coefficients  for  use  by  the  interpolating 
subroutine  SPLINE.  The  arguments  are 

N Number  of  data  (<_100) 

Z Array  containing  values  of  independent  variable 

C Array  containing  values  of  dependent  variable 


\ 

I 


A 


A 


Array  containing  computed  spline  coefficients 
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r 


t 


I 


i 


i 


t 


k 


The  spline  coefficients  are  equal  to  one-sixth  the  second  derivative 
and  are  calculated  by  an  efficient  method  due  to  L.  Solomon  of  Planning 
Systems,  Inc.  A(l)  and  A(N)  are  assumed  to  be  zero. 


h.  SPLINE 

Spline  provides  cubic  interpolation  on  the  basis  of  tabulated  values  of 
a function  and  accompanying  spline  coefficients.  The  arguments  are 


Z 

F 

FZ 


ZM 


N 


I 


Interpolation  point 
Interpolated  function  value 
Interpolated  derivative  value 
Arrays  containing  values  of  the 
independent  variable, 

dependent  variable,  and  spline  coefficients 
Number  of  points 
Pointer  index 


The  pointer  index  is  used  as  in  FINT.  The  calculation 
is  ordinarily  bypassed  to  save  tirm  it  is  implemented 


of  the  derivative 
via  ENTRY  SPLINE2. 


-:>! 
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l . GRAPH 

This  printer-plot  subroutine  depends  on  word  size  and  bit-manipulating 
subroutines  specific  to  the  CDC  6600  and  7600,  and  will  have  to  be  replaced 
for  different  machines.  The  argument  list  is  straightforward  - 

X,  Y Arrays  containing  values  to  be  plotted 

X0,  Y0  Minimum  values,  lower  left 

XSCALE ,YSCALE  Axis  ranges 

N Number  of  points  to  be  plotted 

MODE  See  below 

SYMBOL  Hollerith  character  used  for  plot 

LABEL  Ten-character  hollerith  label 

MODE  = 1 Plots  and  prints  the  function  Y(X) 

MODE  =2  Is  used  for  accumulating  functions  when  plots 

of  more  than  one  function  are  desired 

MODE  =3  Is  used  for  the  last  accumulated  function  to 

print  the  entire  set 

MODE  = 0 Plots  and  labels  the  grid  only,  and  should  be  used 

prior  to  the  first  MODE  = 2 call 


uuuu 
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progpam  7moop ( input, outputs 

COMPUTES  THE  NOPmal-MOHE  I Si  T E W k>  A U • w a V £ OSCILLATIONS  AND 
DISPERSION  RELATIONS  OF  THE  UPPER-OCEAN  THERMOCLINF, 

REAL  A ( 8 ) » NQUN f 8 ) 

DATA  NOUN/6HTCL  P'E  ,8HZMonES,  aHIJTTL/ 

5 prTnt  qq 

qq  FORMAT  (lHl/////30X,qqH***  BOA  PROGRAM  Z MODE  # VERSION  JULY  1 P 7 3 **>* 
*/////) 

1 READ  iO0»(A(!).I»1.8) 
ioo  format ( s a i o ) 

c 

DO  2 I * 1 # 3 

I F ( A ( 2 ) , FO « NOUN ( I ) ) GO  TO  f|0,?0#3ft)#I 

2 CONTINUE 

IF(A(nfEQ.2MGO)  GO  TQ  J 
IF(A(n,EU,«HSTOP)  STOP 
PRINT  98, (A(T), 1st ,8) 

98  FORMAT (32X , 2M*  ,*A101 
GO  TCI  1 
C 

10  CALL  TCLINF(A) 

GO  TO  1 
C 

20  CALL  Z MODES ( A ) 

GO  TO  1 
C 

30  CALL  UTIL(A) 

GO  TO  t 
END 


r»or»  nn  n n n n 


SUBROUTINE  TCLINE(A) 

COMMnN/ETGEN/N,DZ,OZO,7A,ZR,9Nfa00),F(200)  ,F7(200),G(200),0Z(200), 
, HN,ERR 

DIMENSION  A(A),vERR(5)iZZ(UQO),TT(uOO),TAC2oO),TDATA(2nO), 

* ZOATAf200),CYM(fl00) 

DATA  vERR, PI, GRAV.CT1 ,CT?/SHTNPUT,THSET,  WHEXEC,5HPRINT,5HGRAPH, 

* 3, 101592655,9,8?, *,6E*u, , 1 5F*U/ 

DO  1 191,5 

TFCAflJ.EO.VERBd))  GO  TO  f 1 0 , 20 # 30 , 00  # 50 ) , I 
1 CONTINUE 

PRINT  99,  A ( 1 ) 

99  FORMAT(20X,31HltLEGAL  INSTRUCTION,  TCllNF  * ,Aiftf2H  *) 

RETURN 


INPUT  TEMPERATURE  PROFILE 

10  DECODE(10,100,A(U) ) NOATA 

100  FORMAT ( 110) 

READ  101,(ZDATA(!),TDATAfI),Tal,NDATA) 

101  FORMAT { 2F  l ft , 0 ) 
return 

SET  THERMOCLINF  AND  INTEGRATION  PARAMETERS 

20  DECOOE(3O,20O,A(9) ) 7A,ZR,N 
200  FORMAT(2Fl01ft,Ilft) 

oz>za/(n-d 

DZD»nZ**2 

N2«2*N-1 

RETURN 

COMPUTE  INTERPOLATED  STAPJl ITY  PROFILE 

30  CALL  SPLTNFSfNDATA.ZDATA.TDATA.TA) 

J«1 

DO  35  T * 1 ,N2 

Z«DZ*,5*(  I-l) 

7»AMAX1 (Z,7DATA(1 ) 1 

CALL  SPLlNEZfZ.T.OT, 70ATA,T0ATA,  TA.NDATA, J) 
IZ(I)«-Z 
T T ( I ) * T 

ONf N?A1 *I)b«GRAV* (C  T 1 ACT? *T ) *DT 

35  CONTINUE 

SUM9fl  , 

DO  36  1*1, N2 

ON !■ SORT  f AMAXJ  f QN( I ) , 0 , ) ) 
CYN(N2+i-n«ONl#18ftO./PI 

36  SUMaSUM+ONl 
MNa,5*DZ*SUM 
RETURN 


PRINT  TEMPERATURE  AnO  STABILITY  PROFILES 

UO  DECODE  f 1 0 # 1 0 0 * A ( /J ) ) T^KIP 
IFdSKIP.EQ.O)  ISKTP*a 
KOUNTa*! 

DO  «5  !■! »N2, ISKTP 
KOUNTaKOUNT+l 

IF(HnD(KOUNT,50) .EO.OJ  print  uoi 
iiOl  FORMAT(lHi//30<.?7M***  T HE  RMOCL I NE  PROFILE  ***//// 
* 25X,8HDEPTH,  m, ax , 7HTFmd,  C # X , \ l Hn ( Z ) « CY/MR/J 
PRINT  «02»ZZfI)»TTm,CYN(T) 

«0?  FOPMATU7X.3F15,?) 

AI 5 CONTINUE 
RETURN 


GRAPH  PROFILES 

50  DECODE(30*500,A(U) 5 T1,T?,WSCALE 
500  FORMA  T ( JF l 0 , 0 ) 

TSCALE*TI*T2 
MD«I 
SYMalHT 
LBL«7HPRDF I LE 

CALL  GRAPH(TTiZZ.T2,*ZA<TSCALE,ZA,M2,MD,SYM,LBU 
3 Ym*?HN 

CALL  GRAPHfCYN.ZZ,Ot ,-7A;w3CALE, ZA,N?,mD,SYm,lBL) 

RETURN 

END 


non 
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SuARHUT I N£  ZM00tfl(*) 

niMtNSIOM  a(A)#verr(F),x(2ooi,y(2oo) 

COHMnN/FL0/TK(«sn  , TGAHMA  rsi  #?0)  ,TOr.AMHA(Sl,20)  ,FT  (?1,20#2) , 

, FTZ(51i20i2)'Z1.Z?»nk 

COMMON/£IGEN/N,OZ#nzO#  7 A,  ZP.ON(U00I,F{20ft)»F7(20O)#Gf2ft0),f5Z(2G0), 
, HN,FRR 

REAL  K # km a y , norm , KC  Y 

DATA  VERB,FUM,PI/5HIMPIJT,  5HSFT,iJHEYEC,5HPRlNT,5HGRAPH,RHFuNCTJON3» 
, 3.Hi5R2653^/ 

C 

no  1 I a 1 * 5 

TF(A(l),FQtVER«(T))  GO  TO  f 1 0 , 20 , 3ft , «0 . 54 ) , I 
1 CONTINUE 

PRINT  1 0 ft » A ( 1 ) 

100  FQPMAT  ( 1 0 X , 30HILLEGAI  INSTRUCTION,  ZMUOES  **,A10.2H**) 

RETURN 

C 

C 

10  CONTINUE 
RETURN 

SET  WAVENUMBFR  SCALE  ANO  RFSQIU T TON 

20  DEC OO E(50#20ft, *(<!)!  kmax,N*,PRR 
200  FORMAKF10, 0,110,  E10. 0) 

I F ( ERR , EO  » 0 » 1 FRRal.F-S 
KHAX«l002*PI*KMAy 
nK«KMAy/fNK-l ) 

00  21  I K * 1 , NK 

21  TK { IK  > «0K* ( I K« l ) 

return 


non 


■ 
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haIn  COMPUTATION  up  MOnE  FUNCTIONS  and  eigenvalues 

30  DECODE  (40  , 300  * A f <j ) ) mi,m2,Z1,?2 

300  FORMAT (2  1 1 0 # ?F  1 0 ,01 
NF1*N*,5«Z1 /HZ 
NF?«Ni,5*Z?/nz 
PRINT  30 1 , NK 

301  FOPHAT  (1M1  /////2RX#  33H***  FIGENVL  COMPUTATION  TIMES  ***.// 

* 30X,1?H*  RESOLUTION, I 3, 1 «h  valUPS/mOdE  *//// 

* 31X,2&HMQnE  Ayr.  ITFWATIQNS  TImE,/) 

no  3?  MOOE*M( , M2 
9TlMEsSEC0N0(X) 

GAMMA3(MPr)E»t51*PI/HN 
OGaOOG*0, 

AVIT»0, 

no  34  TK»1,NK 
K«IK(IK) 

gamma«gamma*dk*  (dg*  ,5*000 ) 

GAMMAO*GAMMA 

nciaDG 

CALL  EIGFNVL (GAMMA, K,  NORM, ONORM, MODE, MPRP, lEPR,  T T ) 
IFfMERR.EO.O.ANO.IERR.FG.OI  GO  TO  3J 
KCV«500  ,*K/PT 

PRINT  302, MODE. KCY.MFRP, TEPR 

302  FORMAT (//«flX,2RHEIGENVL  FATLFD  TO  CONVERGE  AT,/ 

* 4 5 X , 5HMOOF  , 12, /4SX, I 1 HWAVENUM0FP  ,F  A,?/ 

* 49X,7HMPRR  a , I2/45X,7HIEPR  ■ ,T21 
STOP 

C 

33  p«nqrm/onorm 
DG»K*R/GAMMA 
DDG«OG»DG  J 

IF(IKtEO,l)  nDG»DK*R/GAMM a 

TGAMMAUK,MO,-)E)8rfAMMA 

TDGAMMA(lK,MnDf )aDG 

RNORMbSGRT (ONDBMI 

FI ( IK, MODE, H »F (NF 1 ) /RNQPM 

F I ( IK, MODE, 2 ) sF ( WF2) /RNQOM 

FIZdKiMODF,  1 )bF7(MF1  )/(DZ*RNORm) 

FIZCIK, MODE, 2)»F7(NF?)/(DZ*R NORM) 

AVIT«AVIT^!T 

34  CONTINUE 
C 

TIME * SEC nND(x)«s time 
*VIT»AVIT/NK 

PRINT  303, MODE , t V I T , TIME 

303  FORMAT(32X,  I2,F1?,2,F11  ,<j) 

35  CONTINUE 
RETURN 
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PP I NT  DISPERSION  RELATIONS 

40  DECODE(20,300,A(4))MI ,M2 
no  45  M0D£bM!,m2 
°RINT  400, MODE 

400  PORMATf 1H1//80X,OH***  MODE  ,T2,2SH  DISPERSION  RELATIONS  ***//// 
A 20X,8HK,  CY/KM#4X,0HW,  1/SEC  ,4X,AHW,  CY/HR , 4X , BHC , M/8EC#3X, 

a RHCfi,  M/SEC./) 

00  44  !K»1,NK 
KbTK(IK) 

Cal ,/TGAMMA(lK,MnOF) 

WsCaK 

CG*Ca(1 ,-CAK*TnGAMMA(IK,MODE)  ) 

KCYaSOO , *K/PT 
WHRaJ , 6AC aKCY 
PRINT  401,  KCY,w,wMR#C,CG 

401  F0°MATCUX,5Ei?,5) 

44  CONTINUE 

4?  CONTINUE 
RETURN 

GRAPH  DISPERSION  CURVES 

50  D£CODE(10|500, A(H))  TYPE 

500  FORMATfAlO) 

TF(TYPE.EQ.FUN)  GO  TO  60 
LBLalOHDISPERSION 

DECODE (40,501 ,*(4))  Ml , M2, w SCALE. SK ALE 

501  FORMAT (2110, 2F 10,0) 

SCALE*, 002aPI*SKALE 
MD«0 

CALL  GRAPHIX,  Y,0,,0,,S»<ALE,WSCALE,?0  0,MD,SYM,LBL) 

MD«2 

OKK«KMAX/lR<?; 


non 
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r 


00  55  mqd£*h 1 , M2 
MO0E1 aMon(MOOE, 10) 

ENCODE ( 1 , 502, SYM)  MQ0E1 
502  FORMAT  C I 1 ) 

IF ( MODE  t EG , M2 ) MDa  3 
!K*1 

00  52  Jr  1 ,200 
KaDKK* ( J-l ) 

X ( J ) ■ K 

0AMMA*FINTfK,TK,Tr»AMMA(l,M00£)  , T DG  A MM  A ( \ , MODE  ) , NK  , I K ) 
52  V(J)«1«00,*K/(GAMMA*PP 

CALL  GRAPH(X, Y, 0. , 0 , .SCALE,  wSCALE,2  0 0,MO,SYM,LBL) 

55  CONTINUE 
RETURN 

GRAPH  NORMAL  MODE  FUNCTIONS 

60  OECODE(50,60fl,A(ii))  M j , M? , wui.  T » F SC  A L F .7SCALE 

600  FORMAT(3I10,2F10t0) 

00  65  MOOErMl ,M2 
SCALE* 1 , 

FO*0. 

70«-ZSCALE 
fl  y M ■ l h # 

FNCOOEC 10,601 ,IBL)  MODE 

601  FORMAT ( 5HM0UF  ,1?) 

OELK*KMAX/fMULT-l ) 

00  61  Jol ,N 

61  Y(J)*OZ*(J-l)«ZA 

M0*0 

CALL  GRAPHS,  Y,FO,  70.  SCALE  i Z. SC  ALE,  N,  mo,  SVMiLHU 
M0*2 

1 K ■ 1 

00  65  1*1, MULT 
K-OELK*  (I»1)*,9R) 

-3AMMA»FINTfK,'TK,TGAMMA(l  .MODE)  , T OG  A MM  A ( 1 , MODE  ) , NK  , I K ) 
CALL  EIGENVL (GAMMA #K, norm, nNHRM, mope , mFRR,IERR.IT) 

SC  ALE *FSC ALE* SORT (NORM/ 2 A ) 
F0*«8CALE*rt25*,5*(I-l)/(MiJLT-1  ) ) 

IFd.EO.MULT)  m o r j 

CALL  GRAPMfF, Y.Fn, 70, SCALE,  ZSC ALF.N, mq, SYM.L9L) 

65  CONTINUE 
RETURN 
END 


noon  o non 
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SUBRHUT I NE  EIGFNVL(GAMMA,Kj  NORM,  ("INFIRM , MnDE#MERR#IERR»lT) 

CQHMnN/eiGEM/M#D7#nzG, ZA,2B#GNfafl0l # P (20^) *F7 (200) ,0(200), GZ(200)» 

, HN,fRR 
REAL  K,  NORM 

MgRRaO 

lERHaO 

DGMAVS.6/HN 
SGNsl  , 

TF(MOD(MOOE,2) .EQ,0)  SGNb.1, 

ITERATE  TO  EIGENVALUE 

no  in  1*1,10 

CALL  £IGFNPN(gamma,k,  «,*G) 

I E ( SGN*wG,LE , 0 , ) GO  TO  11 
nGAMMAaW/(2,*GAV|MA*wG) 

IE (AB9(DGAMMA) #GT,0GmAX)  DGAMHAaSIG  N(0GHAX#DGAMHA) 

GAMHA»GAMHA*nGAHM A 

IE  ( AflS  f DGAMHA/GAMMA)  ,(.£  ,FRR)  GO  TO  1? 

10  CONTINUE 

11  TERR*1 

WRONSKIAN  SLOPE  HAS  WRONG  SIGN 
CHECK  FOR  CORRECT  MOPE  NllMRER 

12  M»1 

T T®  I 

EPS«OZQ*W/WG 
PZ(N)aEZ(N1«EPS*GZ(Nl 
00  20  !®2,n 
r c i ) *f  f n*FP.s*G(i) 

FZ  C X ) *F2  f I)-EPfl*GZf n 
20  CONTINUE 
N 1 «N«  1 

00  25  I *2 # N l 

IF((F(I«n*Efn,LT.O.).OR,E(T).EO,0,l  C 

25  CONTINUE 

m£PR  aM-MOPE 

te(mepr,ne,o.or,ierr.eg.i : return 


-48- 


COMPUTE  normal tz at  tom  factor 


0N0RM«*WC*FZ ( N5/OZ 
3UM*0  • 

00  30  T«2,N 
30  3L)M«8UM4F  C I)  **2 

NORMbDZ* ( SUM* ,5*F( 1 )**2+F( t )*D7/fc. 1 
I F ( Z B i EQ , Z A ) RETURN 

TFfK.EQ.O,)  CO  TO  SO 
XaK*(ZB*ZA) 

IF(X,CTt10fl.)  GO  TO  (10 
ExaEXPfX) 

3INH0«(EX»1  ,/F.X)**2*,25 
EXa£X**2 

NQRM* NORM*  (ii?5*(EX-ll/FX5/K»,5*(7B-ZA)5*F(l)**2/SINMQ 

RETURN 

<jo  norm*norm>,s*fut**2/k 


return 

50  NORM*NORMAF(1)**2*fZB«ZA1/3, 

return 

END 


o o r» 
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SUflROUTIME  ETGPNJFMfG4MMA>K|WlWf;) 

CO*MON/E I GEn/N.DZ.OZO.Za.Z*. UN f «00) (F (200) »FZ( 200), G<200) #8Z (200) # 

i HN|FHR 
PEAL  k,ko 

DATA  Cb,C2U/ t\bbbbbbhbbbhbhi  . 0 <1 1 bbb  bbbbbbbb/ 

PQOO»GAMMA(J*X.KQ 

GAMHAQb07U*GAHHA**? 

KQ»DZQ*K**2 

c 

TFfK.EQ.O,}  GO  TO  2 

EXc£XP(,2,ak*(ZB*ZA)  1 

tanhxi  ,-exj  n\ ,*ex) 

P ( 1 ) a T A N H / X 
GU  TO  3 
C 

2 F(l)aZB«ZA 

3 pzm»nz 
R(1)»GZ(lJaO, 

FZZla«F(  \ ) *PO(QNf  1 ) ) 

GZZl ««F ( 1 ) *QN  ( 1 ) 

"0  100  I«2,N 

FlBFd-n 
FZl«FZ(I-n 

o i »g  r : - 1 j 
GZ1*GZ(l-n 

P2»F1*,5*FZl  A,  12*5*FZZ! 

G2«Gt*t5*G7l*t  12«>*GZZ1 
PQ2iPQ(UN( J j ) 

F2»F2«C2<U(FZZ1*P2*PQ2) 

G2«G2-C2<1*  (G7Z1  ♦G2*PQ2*P2*GNf  J)} 

PZZ2»-F2*P02 
GZZ2«-G2*PQ2-P?*GN( J) 

PJ»FUFZUC6*(FZ7U2.*FZ721 
G3«G1*GZUC6*(GZZl*2,*GZZ2) 

PWlaPQ  ( Q Nl  ( J+  1 ) ) 

FZZ3«-F3#PG3 
GZZ3»»G3*PG3«F3*GN( J*l ) 

P(n«F3 
G ( I) aG3 

FZ(I)*FZ1tC6*(PZ71*«.*FZ72»FZZ3) 

GZn)*GZlAC6*(GZ?l*<J,*GZ72*G7Z3) 

FZZIbFZZS 
GZZ1»GZZ3 
too  CONTINUE 


Wa*F f NT 

WG»«G(N)*0ZQ 

RETURN 

FNO 


•A 

* 


I 


I ™ 


; 
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FUNCT10N  FINT(Z,7T,FT,FT7,N, T) 

DIMENSION  ZTf  1),PTU).FT7(1  ) 

5 IFCZ.LT.ZTfin  GO  TO  1 
<1  !F(Z,GT,7T(I*n)  GH  TO  2 

DZ^ZT ( I ♦ 1 ) *ZT ( I ) 

U»(Z«ZT(I) )/nz 

V*1 . -u 

FINT«U*FT(I*| )*V*FT(T)-U*V*( (U«V  ) * (F  T ( T ) -FT  ( !♦  1 ) WDZ*  <U*FTZ  ( !♦  1) 
* -V*FTZ(I))) 

RETURN 

2 T * I ♦ 1 

IFd.lT.N)  GO  TO  R 

I*-l 

RETURN 

t I«I-1 

TFd.GT.t)  GO  TU  3 

I«-l 

RETURN 

END 


SUBROUTINE  SPL I NE S ( N , Z # C , A ) 

SUBROUTINE  SPLINES  CALCULATES  THE  CUBIC  SPLINE  COEFFICIENTS, 
CALLING  ARGUMENTS 

N NUMBER  OF  POINTS,  NOT  TC  EXCEED  100 
Z VECTOR  OF  INPUT  DATA  FOR  INDEPENDENT  VARIABLE 
C VECTOR  OF  INPUT  DATA  FOR  DEPENDENT  VARIABLE 
A COEFFICIENTS  OF  SPLINE 

DIMENSION  Z(i),CCi),A(l),VC100),DC(100),DZ(100) 

A ( l ) = 0 , 

A(N)=0, 

N 1 «N*»  1 


SET  UP  DIFFERENCES, 

DO  12  I a t , N 1 

OZ(I)aZ(lTl)-Z(I) 

DC(I)«CC(I*i)-C(I))/OZ(I) 

12  CONTINUE 
Y(l)30, 

IF(N1,LT,2)  GO  TO  17 
SOLVE  TRIDIAGONAL  SYSTEM, 

DO  15  1=2, Nt 

piv=2,o*(oz(i-i)+oz(in-Dzu-n*Y(i-n 

Y(I)«DZ(I)/P1V 

15  A(I)*(DC(n-DC(I»i)-DZ(I-i)*A(I-l))/PiV 
DO  16  I B*2 , N 1 

IsN+l-IB 

16  A(Z)aA(X).Y(X)«A(lH) 

17  RETURN 


SUBROUTINE  8PUTNF(Z,P*FZ,ZT,FT, A,N, n 
DIMENSION  ZT(l).FT(l)|An) 

MaO 

3 IFfZ.LT.ZTf I))  GO  TQ  1 
U JF(Z.GT,ZTfI*l n GO  TO  2 

DZaZT(Ifl)-ZTd) 

U=(Z-ZT(in/DZ 

VbI ,-U 

F*U*FT(l*i)fV*FTm-OZ**?*il*V*(A(n*(l  .♦V)*A(I*|)*(|t*U)) 
IFfM.EO.O)  RETURN 

FZ»(FT ( 1*1 )*FT ( I) J/OZ*OZ*( *(!)*( 1 .-3.*V**2)-A(I+lJ*Cli»J,*U**2)J 
RETURN 

2 IslM 

IF  f XlLT  f N)  GO  TO  a 

I *»  1 
return 

1 Tal-l 

T F ( I , GE , 1 ) GO  TO  3 

Ta-i 

RETURN 

ENTRY  SPLINE? 

Mat 

GO  TO  3 


n n o 
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! 


SUBROUTINE  GRAPHS  x, y, XO, VO i X SC  ALE, v SCALE  ,N, MODE, SYMBnL,LA0PL  I 

LOGICAL  XOUT,YOUT 

DIMENSION  X( 1 ) , Yf n , A{ u ,31)  , XL ABL ( ^ ) # MASK ( i o ) 

DATA  GRID,PRAME,END3,9LANK/ieH  ♦ , i 0 H«- ♦ , 

* 1 0 H ,,10H  / 

IF (MODE ,GT . 1 ) GO  TO  ? 

INITIALIZE  and  plot  OR  io 


1 DO  3 I a l, SI 
DO  3 K«1 , n 
3 A(K, IJsBLANk 
DO  <|  K a j? , 1 1 
A(K, \ )aFPAME 
a A(K,5l)aFRAME 
DO  5 132,50 
A(1 , IJaENDS 
5 A ( 1 i,  n«ENDS 
DO  6 I ■ t 1 * 4 1 , 1 0 
DO  b Ka2,  11 
b A ( K , DaCRID 
DO  7 K«1 , b 

XLABL(K)axn+0,2*vsCALE*(K,n 

YLABL  bYO^O ,?*VSCALE* (K-i ) 
LINEb61-10*K 
50  FORMAT(FP,3,  1*0 
7 ENCODE! 10,50, A ( 1 ,LlNF) ) YL A«l 
MASKf  10)a<>3 
DO  8 K ■ 1 , 9 

fl  MA3K(10-K)stNt77P,A,3HlET(MA3Kfil-K),6) 
IFfMODE.FQ.OIRETIlRN 


o r»  r»  r»  o r» 
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syhbdl  plot 

2 SYMBOL* SHIFT ( SYMBOL ,-5<0  .A. 77  8 
8YH.SYM80L 
00  9 K»1 , 9 

9 SYMaSYMBnu.O. (8HIFT(SYM,*)<a,,N.77B) 
8YM0OL«8YM 
00  10  I ■ 1 # M 

LINE«50,*(1  ,-(V(T)-Yn)/YSCALF)M.5 
KOLUM»100,*(y(I)«XO)/X8CALF+10.5 
Y0UT»LINF,LT. t .0R.LINE.CT.S1 
XOUTbKOLUH.lt , 1 O.OR.KQLUM.RT. 1 1 0 
iF(xouT,nR. yout)  go  to  10 
KWnR0«(KOLUH.n/10 
KHARaKOLUM«10*KWORO 
KwORObKwORDM 

WORO«A(KWORO#LTNF) ,A,,N,M4gK(KHAR) 
A(KWORO,LINEl»WORO.O. (SYM.A.MASKfKHAR)  ) 

10  CONTINUE 

I F ( MQOE , EO , 2 5 RETURN 

PRINT  ACCUHULATEO  GRAPHS 

PRINT  1000,LABEL#((A(K,L)#Kat#in#L»l#Sn 
1000  FORMAT(1H1///15X,A10/(5X, 11 A10I) 

PRINT  2000, (X^ABLCK) ,K«1 

2000  FORMAT(/6X,6(Fl3.«.7xn 

RETURN 

END 


